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^ : Abstract 

We address the issue of bound state in the two-nucleon system in lattice QCD. Our study is 

^ • made in the quenched approximation at the lattice spacing of a = 0.128 fm with a heavy quark 
OO ■ 

mass corresponding to = 0.8 GeV. To distinguish a bound state from an attractive scattering 
state, we investigate the volume dependence of the energy difference between the ground state 

in 

' and the free two-nucleon state by changing the spatial extent of the lattice from 3.1 fm to 12.3 
fm. A finite energy difference left in the infinite spatial volume limit leads us to the conclusion 
that the measured ground states for not only spin triplet but also singlet channels are bounded. 
Furthermore the existence of the bound state is confirmed by investigating the properties of the 
energy for the first excited state obtained by 2x2 diagonalization method. The scattering lengths 
for both channels are evaluated by applying the finite volume formula derived by Liischer to the 
energy of the first excited states. 

PACS numbers: 11.15.Ha, 11.30.Rd, 12.38.Aw, 12.38.-t 12.38.Gc 
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I. INTRODUCTION 



The strong interaction dynamically generates a hierarchical structure: three quarks are 
bound to form a nucleon with an energy of 1 GeV, and nucleons are in turn bound to form 
nuclei with a binding energy of 10 MeV or so per nucleon. This is a multi-scale physics 
that computational physics should explore, and lattice QCD is responsible for explaining 
the nature of nuclei based on first principles. 

Recently, the present authors have made a first attempt to directly construct the Helium- 
3 and Helium-4 nuclei from quarks and gluons in lattice QCD. In order to control statistical 
errors in the Monte Carlo evaluation of the Helium Green's function as well as quark con- 
tractions whose number factorially increases with the nuclear mass number, calculations 
were carried out at a rather heavy degenerate up and down quark mass corresponding to 
= 0.8 GeV in quenched QCD We successfully confirmed the formation of Helium 
nuclei as a bound state. The key was a systematic change of the spatial size of the lattice 
over a sufficiently wide range that allowed a reliable extrapolation to the infinite volume 
limit. After our finding of the Helium nuclei, NPLQCD Collaboration reported an evidence 
of the H di-baryon bound state in Nf = 2 + 1 QCD at = 0.39 GeV investigating the 
volume dependence of the energy shift from twice of the A baryon mass This was 

followed by HALQCD Collaboration which also presented an evidence of the H di-baryon, 
but in degenerate Nf = 3 QCD at = 0.67-1.02 GeV based on analysis with the effective 
potential measured by the two-baryon wave function ^. 

The situation is markedly different for deuteron. This is the simplest nucleus composed 
of two nucleons in the spin triplet channel, and yet evidence based on lattice QCD for bound 
state has never been reported. It is already quite some time ago that a first analysis of the 
two-nucleon system was made in quenched QCD . Much more recently, studies were 
made with a partially-quenched mixed action [7[ and Nj = 2+1 anisotropic Wilson action 



Extraction of t he p otential between two nucleons has been investigated in quenched and 2+1 
fiavor QCD [SHlQl]. All these studies, however, tried to calculate the two-nucleon scattering 
lengths assuming, based primarily on model considerations with nuclear potentials, that 
the deuteron becomes unbound for the heavy quark mass, corresponding to mj^^ 0.3 GeV, 
employed in their simulations. 

It is time to check the validity of this assumption. We need to investigate whether the 
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bound state exists or not in the heavy quark mass region, where studies so far have been 
carried out, using the arsenal of methods solely within lattice QCD. If there is a bound state, 
the ground state energy never yield the scattering length if substituted into the Liischer's 



finite volume formula 
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12| . In such a case, the scattering length should be obtained from 



the energy of the first excited state. 

We carry out two types of calculations at a heavy quark mass corresponding to = 0.8 
GeV in quenched QCD. The first one is a conventional analysis in which we investigate the 
volume dependence of the energy shift for the ground state. Different volume dependence 
is expected for scattering and bound states. In the second one we investigate the energy 
level of the first excited state employing the diagonalization method [l^ to separate the 
first excited state from the ground state near the threshold of 2m7v. If we find the ground 
state slightly below the threshold and the first excited state slightly above it, then such 
a configuration of the two lowest levels is consistent with the ground state being a bound 
state and the first excited state a scattering state with almost zero relative momentum. 
This method was previously used in a scalar QED simulation to distinguish a system with 
or without a bound state 14\- 

Hereafter we call the analyses employed in the first and second calculations the single 
state and two state analyses, respectively. We also refer to the configuration sets used in the 
two calculations as the first and second ensembles. We should note that the ^Si-^Di mixing 
is neglected in this paper, since we restrict ourselves to measure states in the small relative 
momentum region. 

This paper is organized as follows. Section [III presents the results of the single state 
analysis in the first calculation together with the simulation details. In Sec. IIIII we explain 
the operators employed in the diagonalization method and examine the results obtained by 
the two state analysis. Conclusions and discussions are summarized in Sec. IIVI 



II. SINGLE STATE ANALYSIS 

Let us first present the results of the single state analysis for the ^Si and ^Sq channels. 
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A. Simulation details 

The first ensemble is exactly the same as in the previous work of Ref. We explain the 
parameters once again for clarity. 

We generate quenched configurations with the Iwasaki gauge action [isl at /3 = 2.416 
whose lattice spacing is a = 0.128 fm, corresponding to = 1.541 GeV, determined 



with ro = 0.49 fm as an in put 16| . We employ the HMC algorithm with the Omelyan- 



Mryglod-Folk integrator 



17 



18| . The step size is chosen to yield reasonable acceptance rate 
presented in Table H We take three lattice sizes, xT = 24^ x 64, 48^ x 48 and 96^ x 48, to 
investigate the spatial volume dependence of the energy difference between the two-nucleon 
ground state and twice the nucleon mass. The physical spatial extents are 3.1, 6.1 and 12.3 
fm, respectively. 



We use the tadpole improved Wilson action with csw = 1-378 [16[. Since it becomes 
harder to obtain a reasonable signal-to-noise ratio at lighter quark masses for the multi- 
nucleon system, we employ a heavy quark mass at k = 0.13482 which gives m^r = 0.8 GeV 
for the pion mass and rriN = 1.6 GeV for the nucleon mass. Statistics is increased by 
repeating the measurement of the correlation functions with the source points in different 
time slices on each configuration. The number of configurations and measurements on each 
configuration are listed in Table HI We separate successive measurements by 100 trajectories 
with r = 1 for the trajectory length. The errors are estimated by jackknife analysis choosing 
200 trajectories for the bin size. 

The quark propagators are solved with the periodic boundary condition in all the spatial 
and temporal directions using the exponentially smeared source 

g'(f,t) = X;^e-^''"'\(y,t) (1) 

y 

after the Coulomb gauge fixing. On each volume we employ two sets of smearing parameters: 
{A, B) = (0.5, 0.5), (0.5, 0.1) for L = 24 and (0.5, 0.5), (1.0, 0.4) for L = 48 and 96. The onset 
of ground state can be confirmed by consistency of effective masses with different sources 
as shown later. Hereafter the nucleon operators using the first and the second smearing 
parameter sets are referred to as Oi and O2, respectively. 
The interpolating operator for the proton is defined as 

Pa = £abc{[UaTC'y5dh)u", (2) 



where C = 7472 and a and a, 6, c are the Dirac index and the color indices, respectively. 
The neutron operator Ua is obtained by replacing by rf" in the proton operator. To save 
the computational cost we use the nonrelativistic quark operator, in which the Dirac index 
is restricted to the upper two components. 

The two-nucleon operators for the ^Si and ^Sq channels are given by 

NNssAt) = ^ MM) - , (3) 

NN^s,{t) = ^ Wit)p4t) - p^it)p^it)] . (4) 

For the source operator we insert the smeared quark fields of Eq. ([T]) for each nucleon operator 
located at the same spatial point x. Each nucleon in the sink operator, on the other hand, is 
composed of the point quark fields, and projected to have zero spatial momentum. We call 
this type of sink operator the point sink operator. In the spin triplet channel the operators 
for other two spin components are constructed in a similar way. We increase the statistics 
by averaging over the three spin components. 



B. Numerical results 

Let us first present the effective mass of the nucleon on the (6.1 fm)^ box in Fig. [H 
We observe that the signals with the Oi^2 source operators are clean and the plateaux 
show reasonable consistency with each other. The exponential fit results with one standard 
deviation errors are denoted by the solid lines. They also show the consistency between the 
results from the two nucleon correlation functions. 

Figure [2] shows the effective energy plots for the two-nucleon correlation functions with 
the Oi^2 operators in the ^Si channel on the same volume as in the above. We find clear 
signals up to t ^ 12, beyond which statistical fiuctuation dominates. The effective masses 
with the different sources show a reasonable agreement in the plateau region. The result of 
exponential fit over the plateau region is presented by the solid lines for each operator in 
the figure. Similar behavior of the effective energy is observed in the ^Sq channel as shown 
in Fig. 131 

In order to determine the energy shift AEl = Enn — 2i?Ar precisely in each volume, we 
define the ratio of the two-nucleon correlation function divided by the nucleon correlation 
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m = (5) 



function squared, 

GNN(t) 

{GN{t)y 

where the same source operator is chosen for Gnn(^) and GN(t). The effective energy shift 
is extracted as 

AEf = ln(-^). (6) 



R{t + 1] 

In Fig. m we present typical results of time dependence of AE'j^ for the Oi^2 sources in 
the ^Si channel, both of which show negative values beyond the error bars in the plateau 
region of t = 8-11. Note that this plateau region is reasonably consistent with that for the 
effective mass of the two-nucleon correlation functions in Fig. [2l The signals of AE'j^ are 
lost beyond t 12 because of the large fluctuations in the two-nucleon correlation functions. 
We determine AEl by an exponential fit of the ratio in the plateau region, t = 8-12 for 
Oi and t = 7-12 for O2, respectively. The systematic error of the fit is estimated from the 
difference of the central values of the fit results with the minimum or maximum time slice 
changed by ±1. We obtain a similar quality for the signal for different boxes of (3.1 fm)^ 
and (12.3 fm)'^ as shown in Figs. |5]and[6l respectively. 

The result for the ^Sq channel on the (6.1 fm)^ box is shown in Fig. [71 We find that the 
effective energy shift AE'j^ is negative beyond the error bars, though its absolute value is 
smaller than the ^Si case. The energy shift AEl is determined in the same way as for the 
^Si channel. 

The volume dependence of the energy shift AEl for the ^Si channel is plotted as a 
function of in Fig. |8l Table [TTl summarizes the numerical values of AEl on three 

spatial volumes, where the statistical and systematic errors are presented in the first and 
second parentheses, respectively. The results for the Oi^2 sources are consistent within the 
error bars. Little volume dependence for AEl indicates a bound state, rather than the 
dependence expected for a scattering state, for the ground state in the ^Si channel. 

The binding energy in the infinite spatial volume limit in Table Hi] is extracted by a 
simultaneous fit of the data for the Oi^2 sources employing the fit function including a finite 
volume effect for the two-particle bound state [l^ | 




AEL = --^{i + ^y: -^^ , (7) 
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where 7 and are free parameters, n is three-dimensional integer vector, and denotes 
the summation without \n\ = 0. The binding energy, —AEao, is determined from 7, 

7^ 

-AE^ = -^, (8) 

where we assume 

7^ 



2Jm2,-72_2m^^-^. (9) 

The systematic error is estimated from the difference of the central values of the fit results 
choosing different fit ranges in the determination of AEl, and also using a constant fit as an 
alternative fit form. Adding the statistical and systematic errors by quadrature, we obtain 
— A£'oo=9.1(1.3) MeV for the binding energy. 

We conclude that the ground state in the ^Si channel is a bound state. The provisos, of 
course, are that pion mass is quite heavy and that quark vacuum polarizations are left out. 
Whether these are the reasons why the binding energy is about four times larger than the 
experimental value, 2.22 MeV, is an interesting issue for future study with lighter pion mass 
in full QCD. 

Figure [9] plots the volume dependence of the energy shift AEl for the ^Sq channel, whose 
numerical values are summarized in Table HIl Employing the same analysis as in the ^Si 
channel, we find that —AE^o = 5.5(1.5) MeV in the infinite volume limit, which is 3.7 a 
away from zero. This tells us that the ground state in the ^Sq channel is also bound at 
= 0.8 MeV. Since the existence of the bound state in this channel is not expected at 
the physical quark mass, it might be a consequence of much heavier quark mass used in our 
calculation. Although there are several model calculations varying the up and down quark 



masses, they are restricted around the physical values |20l-|23| . It is an intriguing subject 
to check if the bound state in the ^Sq channel disappears at lighter quark masses. This is 
beyond the scope of this paper, however. 

III. TWO-STATE ANALYSIS 



In this section we present the results of analysis with the diagonalization method |13 |. 
The focus of the analysis is the characteristic feature, well known from quantum mechanics, 
that the existence of a bound state implies a scattering state just above the two particle 
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threshold, and hence a negative scattering length. Our investigation is carried out with the 
diagonalization of 2x2 correlation function matrix. 

A. Simulation details 

We work with two spatial extents, 4.1 fm and 6.1 fm. The corresponding lattice sizes are 
L^xT = 32'^ x48 and 48^ x48, respectively. The latter is the same size as in the first ensemble, 
but we regenerate independent configurations employing the same algorithm. Most of the 
simulation parameters, including the gauge and fermion actions, lattice spacing, quark mass, 
are identical to those explained in Sec. Ull However, the number of configurations and the 
separation of trajectories between each measurement, and the number of measurements on 
each configuration are different. These numbers are tabulated in Table IIIII The errors are 
estimated by the jackknife analysis choosing 400 and 200 trajectories for the bin size on the 
(4.1 fm)^ and (6.1 fm)^ boxes, respectively. These bin sizes are sufficiently large to remove 
the autocorrelation. We use the same operators for the nucleons and two-nucleons as in 
Eqs.([2]), (|3]) and (jl]) choosing the nonrelativistic components. 

The diagonalization method for the 2x2 matrix requires two operators each at source 
and sink time slice, which are explained in the following subsections. 

1. Source operators 

We use the two-nucleon operator composed of the Oi nucleon operator explained in Sec.HTl 
as one of the source operators for diagonalization; from the single state analysis, we expect 
that it has good overlap with the ground state of the two-nucleon system. To reduce the 
statistical error as much as possible we carry out more than one hundred measurements on 
each configuration by changing the center of the smearing source in the spatial and temporal 
directions. 

The diagonalization procedure requires another operator which reasonably overlaps to 
the first excited state. If we envisage this to be a scattering state of the two nucleons with 
almost zero relative momentum, then a possible candidate is an operator consisting of two 
nucleons each projected to zero spatial momentum. Constructing such an operator at the 
source time slice can be done using Z{3) noises for the quark fields. It is empirically known. 
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-1 -^ranc 



however, that statistical noise overwhelms signal in the (multi-)nucleon correlation function 
if the noise is spread over the entire spatial volume. 

The large fluctuation can be reduced by restricting noise to a subset of lattice sites at 
a fixed separation of A'^mod in each spatial dimension at the source time slice. This source 
operator, which we shall call the Or source, is defined by 

M ^ , (10) 

-''rand • , -*^\rt 

where A'rand is the number of the noise, and the color and Dirac indices of the quark field 
are omitted for simplicity. We use the smeared quark fields g'(x, t) of Eq. ([T]) after Coulomb 
gauge fixing, whose parameters are the same as in the 0\ source (A, i?) = (0.5, 0.5) to obtain 
faster plateau of the nucleon state. The smeared quark field is located at 

\/' = {f = xo + nAT^od, ix), < L} (11) 

with xq being a reference position, (xo)i < A^mod? and n being three-dimensional integer 
vector. The complex Z{3) random number ^j{x) satisfies {C,j{x))^ = 1 and has the property 

.r ^™ Tt = ^x,ySg,^. (12) 

The parameters A'mod and A'rand for each calculation are summarized in Table IIIII 



2. Sink operators 

We also need two operators on the sink side to carry out diagonalization. Our idea is to 
employ the solution of the Helmholtz equation in three dimensions for the smearing function 
of the two-nucleon sink operator, 

pi(27r/L)n-r' 

W,2{r-) = C,.J2— (13) 

n 

where is a parameter, r is the relative coordinate between two nucleons, and n is three- 
dimensional integer vector. The overall factor Cg2 is determined from the normalization 
condition |M^g2(ri„ax)| = 1. A similar calculation using the solution of the Helmholtz equation 



was previously reported in Ref. [24 1. 
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In the region of |r| closer to the origin, the smearing function should be modified from 
the free form to take into account the two-particle interaction. One way is to calculate 



the two-particle wave function, as has been done for the two-pion system [25|, l26| and the 
two-nucleon systems js-lDl, and use it as input. We take a simpler alternative of modifying 
the smearing function by hand such that it behaves as a smooth constant function around 
the origin rather than a sharp increase or decrease which occurs for the free form. 

The value of is related to the relative momentum of the two-nucleon state as = 
(27r/L)^ ■ g^, where is not an integer in general due to the finite volume effect of the 
two-particle interaction UJ, ll2[. Since we need two smearing functions, we take a pair of 
values of g^, one around zero momentum q"^ ^ and the other around unit of momentum 
g^ ^ 1 being the simplest choices, and make trial runs to find the optimum values of g^. 
Our optimization criteria are that the effective energy of one of the states is close to ground 
state energy obtained in the single state analysis, and that the two smearing functions have 
significantly different couplings to the ground and first excited states. After several trial 
calculations, we choose g^ = 0.184 and 1.3 for the (4.1 fm)^ box, and g^ = 0.1 and 1.1 for 
the (6.1 fm)^ box. Our smearing functions for both volumes are plotted in Figs, [inland [TTl 

We note that we do not use a negative g^ determined from the bound state which cor- 
responds to an exponentially damped smearing function. In the two-nucleon correlation 
function for the Oi source and such a sink operator, we find that higher excited state 
contributions are not suppressed. Hence such an operator is not suitable for the 2x2 
diagonalization of the ground and first excited states we attempt to carry out. 

Let us finally note that we also employ the point sink operator to carry out the single 
state analysis on the second ensemble for a consistency check with the results from the first 
ensemble. 



B. Results for (6.1 fm)^ box 

We first show the results for the ^Si channel on the (6.1 fm)^ box. Let us begin with 
data for the Oi source. Figure [12] shows the effective energies in the ^Si channel for the 
two smearing function sinks, 1^0.1 and Wi,i, and the point sink P. The effective energy for 
the Wq,i smearing function diverges around t = 8, and rises up from below after t = 12 
due to the sign flip of the correlation function. Thus at least two states contribute to the 
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correlation function overlapping to the operator with different signs. On the other hand, 
the result for Wi,i is close to that for the point sink P. Figure [13] is an expanded view on 
this point. 

Let us now look at effective energies for the Or source in Fig. [TH In this case the VFo.i 
result is close to the point sink result, whereas the Wi^i result is lower. 
We diagonalize the following matrix at each t, 

M(t,to) = C(to)"'C(t), (14) 

where to is a reference time and the 2x2 components of the correlation function matrix 
C(t) are given by 

C^At) = G'^Nit) (15) 

with G^^j^it) being the two-nucleon correlation function using the i {i = Oi,Or) source 
operator and the j {j = Wo.i, W^i.i) smearing function for the sink operator. With a choice 
of to = 6 we determine the two eigenvalues \a(t) (a = 0, 1) of M(t, to) at each t and extract 
the energy of each eigenstate a through Xa(t) = exp{—EL ait — to))- 

The effective energies of the eigenstates obtained from the diagonalization are plotted 
in Fig. [151 The energies for the two states are clearly separated in the plateau region. 
The ground state result is reasonably consistent with the result of the single state analysis 
with the Oi source obtained on the first ensemble, which is expressed by the three solid 
lines in the figure. The first excited state is clearly higher than the ground state, but it is 
much lower than the free case with the lowest relative momentum, whose energy is given by 



l^m}^ + (27r/-L)^ denoted by the single solid line in the figure. 

In order to determine the energy shift as in Sec. [Til we define the ratio of the eigenvalue 
obtained from the diagonalization to the nucleon correlation function squared. 



K{t) 

iGN{t)y 



Rc.it) = (16) 



We also define the effective energy shift of the ratio Ra as 

• Rgjt) 
Ra{t+l) 



AEt = ln(J^]. (17) 
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1. Ground state in the ^Si channel 

Figure [16] shows a compilation of all data for the ground state both from the diagonaliza- 
tion analysis as well as from the single state analysis. The solid circles represent the effective 
energy shift of the ground state A^E^^q using the Oi source in the nucleon propagator in the 
denominator of Eq. f|T6|) . The solid squares are the ones using the Or source in the nucleon 
propagator. The diamonds show the energy shift from single state analysis using point sink, 
but based on the second ensemble. Finally the three lines show the estimated ground state 
energy shift from the single state analysis of the Oi source with the point sink from the first 
ensemble. 

We find it gratifying that the diagonalization results for the ground state exhibit clear 
plateaux over a significant time range extending from t = 7. A somewhat higher value of 
the plateau if one takes the Or source in the nucleon propagator can be traced back to a 
systematic shift in the nucleon effective mass itself, see Fig. [T71 so that the difference should 
be regarded as a measure of systematic error. The plateaux are also consistent with the 
result of the single state analysis from the same ensemble (diamonds), which in turn are also 
consistent with that from the first ensemble (solid lines). 

We determine the central value of the energy shift from the exponential fit of the -Ro(^) 
using the nucleon correlation function of the Oi source with the fit range of t = 7-13. The 
systematic error due to estimate of the threshold 2mN is made from the difference between 
the two results with the Oi and Or sources for the nucleon correlators in the denominator 
of the Rq. The systematic error associated with the fit range is estimated by changing 
the maximum or minimum time slice of the fit range by ±1. Table |IV] summarizes the 
numerical values for the energy shift from the diagonalization analysis AE^ q, and the single 
state analysis AE^. The statistical and systematic errors are presented in the first and 
second parentheses, respectively. We employ an asymmetric systematic error for AEl^q to 
properly reflect an upward shift for the Or source relative to the Oi source. 

2. First excited state in the ^Si channel 

Figure fTSl shows the effective energy shift of the first excited state AE^^. Once again, 
we find a long plateau for both Oi and Or whose values are mutually consistent. The very 
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important feature is that the plateaux are definitively above the threshold and significantly 
lower than the value expected from the free two-nucleon state with unit relative momentum 
as presented by the single solid line in the figure. This is consistent with the ground state 
being a bound state. 

The energy shift AEi^i is determined from an exponential fit of the Ri{t) using the Oi 
source nucleon correlator with the fit range of t = 7-13. We use the fit result for the Or 
source nucleon correlator to estimate a systematic error of the energy shift. The numerical 
value is given in Table |V] . 

3. Analysis of the ^Sq channel 

In the ^So channel, the behaviors of the two-nucleon correlation functions are similar to 
those in the ^Si channel, so that we will present only the results after the diagonalization. 
The effective energies of the eigenstates are shown in Fig. Uni The signals are clean and both 
results show clear plateaux. We observe that the ground state energy is consistent with the 
result on the first ensemble of (6.1 fm)^ box denoted by the three solid lines. 

In Figs. [20] and [21] the effective energy shift for the ground and first excited states AE^^ „ 
are respectively shown as well as the result of AE'^ for the ground state calculated on 
the second ensemble. We find features similar to those in the ^Si channel, including long 
plateaux and systematic biases due to the choice of the source operators. The results for 
energy shift for the ground and first excited states are summarized in Tables HV] and IVt 
respectively, where the errors are estimated as in the ^5*1 channel. 

We observe that the absolute value of the energy shift of the ground state is almost half 
of that in the '^Si channel. This is consistent with the observation in the first calculation. 
On the other hand, the energy shift of the first excited state shown in Fig. |2T]is almost twice 
larger than that in the ^Si channel in Fig. [181 This finding is consistent with the property of 
a system which contains a shallow bound state: The scattering length negatively increases 
as the binding energy decreases, diverging when the binding energy vanishes. 

We confirm then that the two-nucleon system in the ^Sq channel at the heavy quark mass 
of = 0.8 GeV has a bound state as in the ^Si channel. 
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C. Results for (4.1 fm)'^ box 



Scattering states have sensitive dependence on the spatial volume whereas bound states 
do not change much once the spatial size is sufficiently large to contain them. We repeated 
the diagonalization analysis on a (4.1 fm)'^ box to examine if such a difference of the two 
types of states can be conffimed for the ground and ffist excited states in our case. 

The effective energies of the two-nucleon correlation functions with the Oi and Or source 
operators for the ^Si channel are plotted in Figs. [22] and [23l respectively. The behavior we 
observe is similar to the case of the (6.1 fm)^ box except that the effective energy with the 
smearing function Wi,3 has a visible slope in the region where the point sink result shows a 
plateau. 

Figure [23] presents the diagonalization results for the two-nucleon effective energy em- 
ploying to = 8 for the reference time in Eq. iHM . For the ground state, it is once again 
gratifying to find a plateau over a sizable range of time, with the value consistent with 
that from the single state analysis. However, the effective energy for the ffist excited state 
exhibits a visible slope, which was not seen in the (6.1 fm)'^ box case. The positive slope in 
the ffist excited state indicates the presence of contaminations from higher excited state in 
the correlation functions. 

While the present results are not as satisfactory as for the (6.1 fm)"^ box case, we find it 
encouraging that the energy shift relative to the two nucleon threshold, which is negative 
for the ground state, is clearly positive for the ffist excited state and is much lower than the 
value expected for relative momentum of 2-7r/L, see Fig. [25] Because of the presence of a 
positive slope, the estimate of the energy shift suffers from a sizable systematic error from 
the choice of the fit range. We estimate it by making three fits over the ranges t = 9-13, 
11-13, or 12-14, and taking the difference from the ffist one which we use as the central 
value. For the ground state we use t = 9-13 as the central fit range, and shift the minimum 
and maximum time by ±1. The systematic error due to the choice of the Oi or Or source 
is also taken into account. The results are summarized in Tables [IV] and [V]for the ground 
and ffist excited states, respectively, on the (4.1 fm)"^ box for both the ^Si and ^Sq channels. 

In Fig. [26] we plot the energy shift for the ffist excited state from the two lattice volumes 
as a function of A roughly linear behavior, with a larger shift on the (4.1 fm)^ 

box compared to a smaller shift on the (6.1 fm)^ box, is consistent with this state being a 
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scattering state. We evaluate the scattering length using Liischer's finite volume formula 
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I2I ]. and list them in Table \V\ where we find reasonable consistency between the two volumes. 



If our finding of a bound state in quenched QCD at heavy quark mass smoothly continues 
to the physical point, then this is the first calculation which explained a negative scattering 
length for the deuteron channel. 



D. Binding energy from the two calculations 

We evaluate the binding energy of the bound state in the ^Si and ^Sq channels using 
the combined results obtained from both the first and second calculations. Figures [27] and 
281 are the same as Figs. [8] and [H respectively, but including the results of the second 
calculations. The new data are reasonably consistent with the previous ones. We apply the 
same extrapolation procedure to the infinite volume limit as in Sec. [Tll From the fits we 
obtain the following binding energy for the two channels: 

_^^^^,^-H^^^^)m^^^<o.'s,, ^^^^ 



4.4(0.6)(1.0) MeV for ^S, 



0) 



where the first and second errors are statistical and systematic. These results are reasonably 
consistent with the ones in Sec. [Tll 



IV. CONCLUSION AND DISCUSSION 



We have carried out two calculations in quenched QCD to investigate whether the two 
nucleon systems are bound or not at the heavier quark mass, corresponding to m,r = 0.8 
GeV. In the first calculation, we have focused on the ground state of the two-nucleon system, 
and have investigated the volume dependence of the energy shifts obtained with two different 
source operators. We have found that the ground state in the ^Si channel has little volume 
dependence, and a finite energy shift remains in the infinite volume limit. Based on these 
results we have concluded that the ground state is a bound state at the heavy quark mass. 
A similar result is obtained in the ^Sq channel, though the binding energy is almost half of 
the one in the ^Si channel. 

In the second calculation we have carried out two-state analysis using the diagonalization 
method. The ground and first excited states are well separated on the (6.1 fm)^ box, and the 
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ground state energies for the two channels agree with the ones obtained from the single state 
analysis. The energy of the first excited state is positive and far below the free two-nucleon 
energy with the lowest relative momentum in both channels. This leads to the conclusion 
that each channel has one bound state. We obtain similar results on the (4.1 fm)'^ box, 
though the contaminations from higher excited states may be larger than the (6.1 fm)^ case. 
The energy of the first excited state increases as the volume diminishes. The scattering 
length is obtained from the energy of the first excited state using the finite volume formula. 
The results in the two volumes reasonably agree with each other. In the ^Si channel the 
scattering length is roughly one fifth of the experimental value. The difference might be 
attributed to the heavier quark mass employed in this calculation. 

The existence of the bound state and the negative scattering length in the ^Sq channel 
looks odd from the experimental point of view. In addition we cannot directly compare 
our result with those of the model calculations, which are restricted around physical quark 
masses. We expect that the bound state vanishes at some lighter quark mass, where the 
scattering length diverges changing the sign from negative to positive. Further reduction of 
the quark mass would decrease the scattering length. Confirmation of this scenario requires 
to investigate the quark mass dependences of the binding energy and the scattering length. 
We leave this study to future work. 
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TABLE I: Number of configurations (A'^conf)) number of measurements on each configuration 
(-^meas), acceptance rate in the HMC algorithm, pion mass (m7r) and nucleon mass {m]\f) for the 
first ensembles. 

L iVconf iVmeas aCCept.(%) 711^ [GcV] niN [GeV] 

24 2500 2 93 0.8000(3) 1.619(2) 

48 400 12 93 0.7999(4) 1.617(2) 

96 200 12 68 0.8002(3) 1.617(2) 



TABLE II: Energy shift —AEl in MeV units for ^Si and ^Sq channels on each spatial volume with 
the first ensembles. Extrapolated results to the infinite spatial volume limit are also presented. 
The first and second errors are statistical and systematic, respectively. 







'Si 




L 


Oi 


O2 


Oi O2 


24 


10.2(2.2)(1.6) 


10.0(1.5)(0.5) 


6.1(2.3)(2.2) 8.4(1.5)(0.5) 


48 


9.6(2.6)(0.9) 


10.2(2.0)(0.8) 


5.2(2.6)(0.8) 6.4(2.0)(0.8) 


96 


7.8(2.1)(0.4) 


9.0(2.0)(0.5) 


4.6(2.0)(1.1) 6.0(1.9)(0.5) 


oo 


9 


.1(1.1)(0.5) 


5.5(1. 1)(1.0) 



TABLE III: Number of configurations {Nconf)^ separation of trajectories between each measure- 
ment (A'sep), number of measurements with Oi source on each configuration (Nj^cas), number of 
Z{3) random number for Or source on each configuration (iVi-and)) spatial interval between smeared 
quark fields for Or source {N^q^), acceptance rate in the HMC algorithm, pion mass (771,^) and 
nucleon mass (ttitv) for the second ensembles. 



L 


^conf 


A^sep 


A'rncas 


Arand 


A'mod 


accept. (%) 


[GeV] 


mpf [GeV] 


32 


300 


400 


192 


40 


16 


87 


0.7998(2) 


1.6162(9) 


48 


300 


200 


144 


32 


12 


93 


0.8001(1) 


1.6176(4) 
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TABLE IV: Energy shifts —AEi and —AEifi in MeV units for '^Si and -"^Sq channels at L = 32 and 
48 on the second ensembles. The first and second errors are statistical and systematic, respectively. 



L 


-AEl 


-AEl,o 


-AEl 




32 

48 


7.9(0.6)(0.8) 
8.5(1. 1)(0.3) 


6.4(1.3) (i°:D 
7.1(0.7) (l°i) 


4.7(0.7)(0.6) 
4.8(1.0)(0.7) 


3.0(1.7)(i°;D 
4.5(0.9) (t°:l) 



TABLE V: Energy shift of the first excited state AEl^i and scattering length ao for ^Si and ^Sq 
channels at L = 32 and 48 after diagonalization in two-state analysis. The first and second errors 
are statistical and systematic, respectively. 













L 


AEl,i [MeV] 


«o [fm] 


AEl,i [MeV] 


ao [fm] 


32 


13.3(1.3) (tf^) 


-1.5(0.2) (t?i) 


15.8(1.6)(t^:^) 


-1.8(0.3) 


48 


2.3(0.8) (1^:?) 


-1.05(24) {tlil) 


4.2(0.8) 


-1.62(24) (t°:?^) 
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FIG. 1: Effective mass of nucleon with Oi (circle) and O2 (square) sources on (6.1 fm)^ box in 
lattice unites. Fit results with one standard deviation error band are expressed by solid lines. 
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FIG. 2: Effective energy of two- nucleon for ^Si channel with Oi (circle) and O2 (square) sources 
on (6.1 fm)"^ box in lattice unites. Fit results with one standard deviation error band are expressed 
by solid lines. 
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FIG. 3: Same as Fig. [2] for ^Sq channel. 
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FIG. 4: Effective energy shift of two-nucleon AE*^ for ^Si channel with Oi (circle) and O2 
(square) sources on (6.1 fm)'^ box in lattice units. Fit results with one standard deviation error 
band are expressed by solid lines. 
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FIG. 5: Same as Fig. Hon (3.1 fm)^ box. 
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FIG. 6: Same as Fig. Hon (12.3 fm)^ box. 
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FIG. 7: Same as Fig. [4] for ^Sq channel. 
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FIG. 8: Spatial volume dependence of AE^ = — 2mj\[ in GeV units for ^Si channel with 

Oi (circle) and O2 (square) sources. Statistical and systematic errors are added in quadrature. 
Extrapolated results to the infinite spatial volume limit (filled circle) and experimental values (star) 
are also presented. 
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FIG. 9: Same as Fig. [8] for ^Sq channel. 
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FIG. 11: Same as Fig.IOlon (6.1 fm)^ box. 



2.5 



2.4 - 



2.3 



2.2 



2.1 - 



2.0 - 



Li 



1.9 



O 

□ w,, 
O P 



<D 



(D ^ o 



8 12 
t 



16 



20 



FIG. 12: Effective energy of two-nucleon for ^Si channel with d source on (6.1 fm)^ box. Results 
are given in lattice units with smearing functions VFq.i (circle) and VFi.i (square), and point sink 
operator (diamond). 
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FIG. 13: Same as Fig. [T2\ but scale in vertical axis is enlarged. Solid lines denote fit result for 
ground state energy with one standard deviation error band in single state analysis choosing Oi 



source. 
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FIG. 14: Same as Fig. [T3l with Or source. 
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FIG. 15: Effective energy of two-nucleon ground (circle) and first excited (square) states obtained 
by diagonalization method for '^Si channel in lattice units. Fit result of single state analysis with 
Oi source on the first ensemble with one standard deviation error band is expressed by three solid 
lines. Expected energy level of free two-nucleon state with lowest relative momentum is denoted 
by single solid line. 
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FIG. 16: Effective energy shifts of two-nucleon A^E^ g for ^Si channel on (6.1 fm)^ box in lattice 
units. Nucleon correlators with Oi (circle) and Or (square) sources are employed in the denom- 
inator of Eq. (jl6p . Result of single state analysis with Oi source AE'j^ is also plotted by open 
diamond symbol. Fit result of single state analysis with Oi source on the first ensemble with one 
standard deviation error band is given by solid lines. 
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FIG. 17: Nucleon effective mass on (6.1 fm)^ box for the second ensemble with Oi (circle) and Or 
(square) sources in lattice units. Fit result of Oi source on the first ensemble with one standard 
deviation error band is expressed by solid lines. 
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FIG. 18: Same as Fig. [16] for the first excited state. Expected energy level of free two-nucleon 
state with lowest relative momentum is denoted by single solid line. 
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FIG. 19: Same as Fig. [15] for ^Sq channel. 
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FIG. 20: Same as Fig. [16] for ^Sq channel. 
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FIG. 21: Same as Fig. [18] for ^Sq channel. 
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FIG. 22: Effective energy of two-nucleon of Oi source on (4.1 fm)^ box for ^Si channel in lattice 
units with smearing function Wi.s (square) and point sink operator (diamond). Fit result of single 
state analysis with Oi source on the first ensemble of (6.1 fm)^ with one standard deviation error 
band is denoted by solid lines. 
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FIG. 23: Same as Fig. [22] with Or source. 
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FIG. 24: Same as Fig. [15] on (4.1 fm)^ box. Result of single state analysis is denoted by diamond 
symbol. 
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FIG. 25: Same as Fig. [H] on (4.1 fm)^ box. 
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FIG. 26: Spatial volume dependence of AE^^i in GeV units for the first excited states in ^ 
(circle) and ^Sq (square) channel. The square symbols are slightly shifted to positive direction 
horizontal axis for clarity. Statistical and systematic errors are added in quadrature. 
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FIG. 27: Spatial volume dependence of AE^ and AEl^q in GeV units for ■^Si channel with 
Oi (circle), O2 (square) sources on the first ensembles and Oi source on the second ensemble 
(diamond). Result of AEj^ q obtained by diagonalization method is denoted by triangle symbol. 
The diamond and triangle symbols at ~ 10~^ are slightly shifted to positive direction in 

horizontal axis for clarity. Statistical and systematic errors are added in quadrature. Extrapolated 
results to the infinite spatial volume limit (filled circle) and experimental values (star) are also 
presented. 
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FIG. 28: Same as Fig. EH for ^Sq channel. 
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